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1.1.  Photoelectron  Energy  Spectrum 


In  the  daytime  midlatitude  lower  ionosphere  the  primary  source  of  electrons 
is  photoionization  of  the  neutral  gas  by  photons  in  the  extreme  ultraviolet 
(EUV)  region  of  the  solar  spectrum.  We  have  studied  this  process  in  detail 
along  with  the  various  particle-particle  collisional  processes  that  determin.* 
the  energy  dependence  of  the  photoelectron  flux.  Our  analysis  is  based  on 

Boltzmann-Fokker-Planck  theory  in  the  local  approximation  as  developed  by 
Jasperse* . 


Using  this  theory  we  can  calculate  the  isotropic  electron  distribution 
function  (EDF)  from  a  detailed  balance  equation  of  electron  sources  and  sinks  in 
phase  space.  The  equation  in  energy  space  is 


0 


SFq(E)  SFq(E) 

6t  pi  +  6t 


coll 


(1) 


where  the  terms  on  the  right-hand  side  are  given  in  Jasperse^*^.  The  two 
general  features  of  the  EDF  are  a  Maxwellian  in  the  thermal  region  winb  a  struc¬ 
tured  and  inflated  tail  at  high  energies  (1-200  eV).  The  electrons  are 
generally  produced  at  high  energies  by  photons  ionizing  the  neutrals .  These 
photoelectrons  then  lose  energy  through  both  electron  conserving  and  electron 
producing  collisions  and  eventually  attain  "thermal”  energies  where  recom¬ 
bination  with  ions  becomes  an  effective  sink. 


In  1.1.1,  we  discuss  the  detailed  processes  that  lead  to  the  high  energy 
tail  of  the  EDF.  In  Section  B,  we  focus  on  the  role  of  the  electron-electron 
collisions  in  the  transition  region  of  near-thermal  energies  where  the  EDF 
changes  character  from  a  tail-like  structure  to  a  Maxwellian. 


1.1.1.  High  Energy  Tall  of  the  Photoelectron  Distribution  Function 

In  the  photoionization  process,  a  neutral  particle  absorbs  a  photon  of 
energy  hv  and  an  electron  is  ejected  with  energy  E*hv  -Epijm,  where  Epijm  is  the 
ionization  energy  associated  with  the  ionization  process  in  which  the  neutral  j 
is  transformed  into  an  ion  in  state  m.  The  number  of  electrons  with  energy  E 
produced  by  such  a  process  is 


I  (E+Epijjj)  Qpi  jm(E+Emi  jm)  Oj  , 


(2) 


where  I  is  the  number  of  photons  with  energy  hv=E+Epijm,  Qpijm  Is  the  cross- 
section  for  the  photoionization  process  described  above,  and  nj  is  the  density 
of  neutral  j.  To  calculate  the  total  electron  production  rate  due  to  a  photon 
of  energy  hv  one  has  to  sum  over  all  the  neutral  species  and  all  the  final  ionic 
states.  In  our  initial  modeling  we  included  contributions  from  20  of  the 
electronic  bound  and  dissociating  states  of  0+,  N2+  and  02+. 

The  solar  EUV  spectrum  is  highly  structured  being  composed  primarily  of 
narrow  lines'^  with  the  largest  line  being  He  II  (40.7  eV)  with  a  measured  width 
(FWHM)  of  .013  eV.  From  formula  2,  it  is  clear  that  if  the  cross  section  is 
reasonably  constant  over  the  width  of  a  solar  line  then  such  a  narrow  line  at  hv 
will  lead  to  a  narrow  peak  in  the  photoelectron  production  spectrum  at  hv-Epijm. 
Such  a  production  peak  is  then  degraded  by  various  interactions  in  the  plasma 
resulting  in  shifted,  reduced  and  smeared  peaks  in  the  photoelectron  spectrum. 
Such  structure  has  been  seen  in  measured  photoelectron  spectrum^,  particularly 
between  20  to  30  eV  where  several  large  peaks  can  be  identified  with  the  photo¬ 
ionization  [of  O2  and  N£]  by  the  He  II  line. 

In  spite  of  the  particle  interactions  much  of  the  structure  In  the  photo¬ 
electron  spectrum  can  be  directly  Identified  with  a  specific  solar  line,  neutral 
specie  and  final  ion  state.  The  four  major  peaks  of  the  photoelectron  spectrum 
form  a  definite  pattern  determined  by  threshold  energies  for  producing  N2+  and 
02+  in  their  ground,  first  excited  and  second  excited  states.  While  the  He  II 
(40.8  eV)  solar  line  produces  the  largest  version  of  the  pattern  between  21  and 
28  eV  strong  solar  lines  at  33.6  eV  and  48.2  eV  produce  the  same  pattern  simply 
shifted  in  the  photoelectron  spectrum  by  the  energy  difference  between  these 
solar  lines  and  the  HE  II  line.  The  pattern  is  not  limited  to  four  peaks  but  is 
part  of  a  larger  pattern  which  Includes  the  photoproduction  of  other  ionic  states 
by  a  particular  solar  line. 

While  we  can  identify  some  of  the  structure  in  the  photoelectron  spectrum 
as  arising  directly  from  photoionization  there  is  another  process  which  contri¬ 
butes  significantly  to  the  structure,  namely,  electron-neutral  inelastic  scat¬ 
tering.  In  inelastic  scattering  the  electrons  of  the  photoionization  peaks 
scatter  off  the  neutrals  and  lose  discrete  amounts  of  energy  resulting  in 
shifted  peaks  which  are  reduced  images  of  the  photoionization  peaks.  With  these 
two  processes,  photoionization  and  electron  neutral  inelastic  scattering,  we 
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found  excellent  agreement  with  satellite  data  with  respect  to  the  locations  of 
peaks  and  the  general  character  of  the  photoelectron  spectrum. 

While  we  included  such  interactions  as  electron-electron  and  electron- 
neutral  to  account  for  the  degradation  of  the  production  peaks  a  variety  of 
effects  which  would  broaden  a  peak  during  the  photoionization  process  were 
neglected.  To  account  for  this  the  solar  lines  were  given  an  increased  width,  2 
eV  for  HE  II,  where  the  width  was  on  the  same  order  as  the  vibrational  spacing 
for  the  molecules.  This  ad  hoc  width  fundamentally  determined  the  width  and 
shape  of  the  peaks  in  the  photoelectron  spectrum.  While  our  agreement  with  data 
was  generally  excellent  the  widths  of  the  theoretical  peaks  were  less  than  the 
observed  widths.  Before  one  can  consider  if  other  processes  such  as  wave-particle 
interactions  could  be  important  it  is  necessary  to  first  consider  the  photoioni¬ 
zation  process  in  greater  detail. 

We  have  found  that  the  effects  in  photoionization  due  to  natural  line 
broadening,  Doppler  broadening,  pressure  broadening,  molecular  rotation  levels 
and  atomic  fine  structure  are  of  the  same  order  or  smaller  than  the  experimental 
widths  of  the  solar  lines.  From  photoelectron  spectroscopy  data**  we  have  seen 
that  enough  of  the  vibrational  levels  of  a  given  molecular  electronic  bound 
state  can  be  populated  so  as  to  cause  a  spread  in  electron  energy  as  large  as 
4.94  eV.  We  include  these  facts  in  our  modeling  by,  first,  using  a  solar 
spectrum  with  line  width  (.04  eV  for  HE  I)  of  the  order  of  the  measured  solar 
line  widths  and,  second,  explicitly  including  the  molecular  vibration  levels  In 
the  photionization  process. 

Inclusion  of  vibrational  effects  changes  Eq.  (2)  in  the  following  manner 

I(hv)  Bn jm(hv)  Qpijm(hv)  nj  (3) 

where  Bnjm's  are  branching  ratios  for  the  vibrational  states  for  an  electronic 
state  jm7 .  The  electron  energy  E  is  now  given  as  E=hv-Epi jm-^Env^b  so  that  the 
number  of  photoelectron  peaks  from  a  particular  solar  line  and  molecular  state 
jm  equals  the  number  of  vibrational  levels  of  state  jm  that  are  populated.  In 
our  current  modeling  we  have  included  branching  ratios  for  8  of  the  02+  and  N2+ 
electronic  states. 

These  modifications  broaden  the  theoretical  widths  and  remove  the  ad-hoc 
nature  of  our  previous  calculations.  The  energy  resolution  of  present  photo¬ 
electron  data  is  such  that  it  is  difficult  to  quantify  what  differences  remain 


between  theoretical  and  experimental  widths.  It  does  appear  though  that  theory 
may  still  underestimate  the  broadening  of  photoelectron  peaks.  Theoretically 
the  next  step  is  to  include  the  interactions  of  electrons  with  the  thermal  wave 
fluctuations.  Experimentally,  high  resolution  measurements  of  the  tail  of  the 
EDF  would  clearly  quantify  differences  between  any  theoretical  modeling  and 
data. 


1.1.2.  Photoelectron  Distribution  Function  at  Thermal  and  Near  Thermal  Energies 

In  recent  years  there  has  been  increased  interest  in  the  thermal  and  near 
thermal  energy  regions  of  the  photoelectron  distribution  function.  In  theoreti¬ 
cal  calculations  of  the  EDF  how  one  treats  the  electron-electron  interactions 
(at  these  energies)  is  an  important  concern.  From  an  examination  of  our  solu¬ 
tions  of  Eq.  (1)  and  the  various  processes  involved  at  low  energies  we  have  come 
to  some  conclusions  concerning  how  the  electron-electron  interactions  should  be 
treated . 

We  used  the  Fokker-Planck  description  for  the  Coulomb  collisions  between 
electrons 
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The  first  term  is  called  the  dynamical  friction  term  and  the  second  is  the 
velocity  diffusion  term.  The  friction  term  describes  the  tendency  of  electron- 
electron  collisions  to  act  as  a  force  accelerating  or  decelerating  electrons  to 
the  average  velocity.  The  diffusion  term  describes  the  tendency  of  the  colli¬ 
sions  to  diffuse  the  electrons  out  in  velocity  space  so  there  are  random  fluc¬ 
tuations  about  the  average  velocity.  A  Maxwellian  distribution  of  electrons 
represents  a  distribution  where  these  two  effects  exactly  balance. 

Our  calculation,  done  in  energy  space  following  the  Rosenbluth,  MacDonald 
and  Judd^  development  and  involving  the  integral  of  the  kinetic  equation,  gives 

47T  -4  6't~  dE  =  “kee  P7T  "  4lTkee  (f  E1/2(l2+J-i)  F0)  (5) 

where 
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T0  =  4it  F(x)  dx  ,  I2  ■  — g-  xF(x)  dx  ,  J_i  =  4x  E1/2  dx 

3/2  2 

and  kee  =  2(y)  4u  (~)  InA .  The  first  term  we  still  label  friction  and  the 
second  diffusion.  Further,  the  average  energy  loss  rate  for  an  electron  is 

which  in  the  case  of  having  F(x)  a  Maxwellian  in  Iq,  I2  and 

agrees  with  the  Butler  and  Buckingham  energy  loss  rate.^  It  is  some  form  of 
friction  term  which  is  used  in  most  calculations  of  the  photoelectron  distribu¬ 
tion  function  while  the  diffusion  term  is  often  neglected. 

The  use  of  a  friction  term  can  be  justified  in  general  for  high  energy 
where  fast  electrons  are  decelerated  towards  thermal  regions.  The  two  excep¬ 
tions  where  ignoring  diffusion  can  be  a  problem  are  (1)  near  any  sharp  peaks 
where  derivatives  of  F0  can  be  large  and  the  proper  description  of  the  spreading 
of  the  peaks  may  require  diffusion  and  2)  in  the  thermal  and  near  thermal 
regions  where  the  F  is  or  is  becoming  Maxwellian  in  character  so  that  the  very 
existence  of  the  Maxwellian  indicates  that  friction  and  diffusion  are  of  equal 
importance. 

[From  our  calculations]  we  have  seen  [in  3  ways]  the  importance  of  the  dif¬ 
fusion  term.  First,  a  comparison  of  our  full  solution  to  a  simpler  CSD  solu¬ 
tion.  The  CSD  solution  Is  an  example  of  a  calculation  where  only  a  friction 
type  term  is  used,  and  we  find  that  while  It  follows  the  high  energy  tail 
without  peaks  that  at  thermal  and  near  thermal  energies  it  contains  none  of  the 
Maxwellian  character  of  the  full  distribution  function.  Second,  we  did  runs 
with  the  diffusion  term  turned  off  above  various  threshold  energies.  When  the 
threshold  was  above  10KT  away  from  the  Maxwellian  the  solution  was  little 
affected.  As  the  threshold  was  lowered  below  10KT  into  the  Maxwellian  region 
the  Maxwellian  character  above  the  threshold  was  completely  lost  and  the 
electron  temperature  became  radically  wrong.  Finally  we  were  able  to  examine 
term  by  term  the  various  processes  included  In  our  calculations.  We  found 
directly  that  at  high  energies  the  friction  term  is  generally  much  larger  than 
the  diffusion  term,  but  as  one  goes  down  to  near  thermal  energies,  10  KTe,  the 
diffusion  terra  becomes  of  the  same  order  as  the  friction  term,  but  opposite  sign. 
Further  the  two  electron  terms  together  do  not  dominate  over  the  other  processes 
of  these  energies,  but  the  friction  term  alone  dominates  the  other  processes  so 
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that,  if  the  diffusion  term  is  not  there  to  cancel  with  it,  a  large  error  would  be 
introduced  in  the  equation. 

Our  basic  conclusion  is  that  the  Maxwellian  forms  a  natural  boundary 
where  one  should  use  both  friction  and  diffusion  to  describe  the  electron- 
electron  interaction.  In  other  words  if  one  is  concerned  with  the  lower 
energies  in  the  transition  region  between  the  Maxwellian  and  the  high  energy 
tail  then  electron-electron  diffusion  effect  should  be  considered  along  with 
electron-electron  friction.  Recent  work*0  has  further  confirmed  the  impact  of 
and  need  for  both  types  of  terms  in  the  calculation  of  the  thermal  electron 
heating  rates. 

1.2.  Plasma  Instabilities  and  Their  Effects  on  the  Photoelectron  Energy 
Distribution 

One  interesting  feature  of  the  theoretically  calculated*  photoelectron 
energy  distribution,  which  is  in  disagreement  with  the  measurements,  occurs  in 
the  2-6  eV  energy  range.  In  this  energy  range  the  distribution  function  has  a 
minimum  at  about  2.3  eV  (Fig.  1-1).  This  is  explained  theoretically,  since  the 
cross  section  for  the  electron-impact  excitation  of  vibrational  states  of  N2  has 
a  maximum  at  2.3  eV.  Beyond  2.3  eV,  the  photoelectron  flux  rises  sharply  to  a 
maximum  at  about  5  eV,  as  the  excitation  cross  section  decreases  sharply. 

Beyond  this  energy  the  flux  falls  off  as  the  cross  section  for  the  excitation  of 
melastable  states  of  atomic  oxygen  increases.  The  theoretical  calculations 
(Fig.  1-1)  show  that  this  minimum  in  the  spectrum  is  more  and  more  pronounced  at 
or  below  altitudes  of  130  km.  Above  130  km,  the  valley  starts  to  be  filled  up, 
the  peak-to-valley  ratio  decreases,  until  at  or  above  210  km  the  structure 
disappears  completely.  This  disappearance  can  be  attributed  to  the  depletion  of 
N2  as  well  as  to  the  gradual  smoothing  process  arising  from  electron-electron 
collisions  as  the  altitude  increases. 

Measurements  of  the  electron  energy  distribution  by  McMahon  and  Heroux**, 
who  studied  specifically  the  2-5  eV  energy  range  with  improved  energy  resolution 
of  the  apparatus,  are  in  good  agreement  with  the  theoretical  calculations  of 
Jasperse*  at  and  above  170  km.  Below  170  km,  the  calculated  values  of  the  peak- 
to-valley  ratios  by  Jasperse*  is  larger  than  the  measured  values.  The  measure¬ 
ments  are  in  most  striking  disagreement  with  the  theory  below  130  km,  where  they 
show  plateaus  in  the  distribution  functions  in  the  2-5  eV  energy  range. 
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This  discrepancy  between  theory  and  measurement  suggests  that  in  the  low 
altitude  regions  (100-170  km)  collisional  processes  alone  cannot  explain  the 
observed  photoelectron  distributions.  It  is  well-known  that  a  homogeneous 
plasma  in  a  magnetic  field  with  isotropic  distribution  functions  can  be 
unstable  if  a  population  of  high  energy  particles  is  also  present.  This  is  pre¬ 
cisely  the  situation  in  the  lower  ionosphere,  and  it  can  be  expected  that 
the  plasma  instability  will  produce  the  anomalous  diffusion  in  the  velocity 
space  through  wave-particle  interaction,  which  in  turn  will  flatten  the  distri¬ 
bution  functions  in  the  2-5  eV  energy  range.  With  this  in  mind,  we  studied 
excitation  of  electrostatic  instabilities  in  the  ionospheric  collisional  plasma 
by  the  suprathermal  electrons  near  the  5-eV  maximum. 

Using  the  Bhatnagar-Gross-Krook  (BGK)  collision  operator  to  take  into 
account  the  effects  of  electron-neutral  collisions,  we  find  that  the  suprather¬ 
mal  electrons  near  5  eV  can  excite  electron  cyclotron  modes  in  the  110-130  km 
region,  and  upper  hybrid  modes  above  125  km^,  As  the  waves  grow  in  amplitude 
at  the  expense  of  the  suprathermal  electron  energy,  growth  rates  are  reduced. 
Eventually,  a  steady-state  is  reached  in  which  the  waves  saturate  to  a  finite 
level  of  amplitude  and  the  electron  energy  distribution  in  the  2-6  eV  energy 
range  is  considerably  modified  through  anomalous  diffusion  in  velocity  space. 

We  have  studied  these  wave-particle  effects  within  the  framework  of  the  so- 
called  quasilinear  theory  of  plasma  turbulence. 

In  this  section*  we  first  discuss  the  linear  theory  of  the  plasma  instabili¬ 
ties  mentioned  above.  We  then  present  the  relevant  quasilinear  equations  and 
discuss  the  approximate  analytic  solutions  of  these  equations,  showing  the  time 
evolutions  of  the  electron  distribution  function  and  the  wave  spectral  energy 
density  from  their  initial  to  the  time-asymptotic  values. 

1.2.1.  Linear  Theory  of  the  Plasma  Instabilities 

The  equation  for  the  one  particle  photoelectron  distribution  function,  f, 
is 

[57  +  I  '  4  -  »  +  '  if] f  '  S  +  1(0  ’  (6) 

where  Bp  is  the  geomagnetic  field  (assumed  to  be  uniform),  S  is  the  photo¬ 
electron  source  function,  and  L  is  the  total  collision  operator  for  electron- 
neutral,  electron-electron  and  electron-ion  collisions.  The  other  quantities 
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have  their  usual  meanings.  We  Imagine  that  the  equilibrium  electron  population 
is  composed  of  a  Maxwellian  part  and  a  suprathermal  part,  denoted  by  the  two 
values  of  the  subscript  j,  and  study  electrostatic  perturbations  from  this 
equilibrium  having  frequencies  u.~u/ p,  where  o»p  is  the  electron  plasma  frequency. 
At  these  frequencies  the  ion  dynamics  can  be  ignored.  We  seek  solutions  of  the 
form 

fj  =  foj(v>  +  flj(l»l»t)  »  (7) 

E  =  0  +  E^r.v.t)  (8) 

where  f0j  is  the  equilibrium  (isotropic)  distribution  of  the  j  type  electrons, 
flj  is  the  corresponding  small  amplitude  perturbation,  and  is  the  perturbed 
electric  field.  Approximately  decoupling  the  thermal  from  the  suprathermal 
populations  we  find  that  the  zero  order  equations  are 

0  -  S  +  L(foj)  ,  (9) 

and  the  first  order  equations  are 

+  l  •  a!"  ~  (--o)  *  41  flJ 
*  J  h  '  4  foj  =  L(f°j+flj>  "  L(f°j)  *  C10> 

g|  *  -j/d3v  flJ  '  (ll) 

It  can  be  shown  that  L  may  be  approximated  by  the  BGK  operator  for  the  dominant 
electron-neutral  collisions,  i.e., 

R.H.S.  of  (10)  =  “v jf 1 j+v j(ni j/n0j)f0j  ,  (12) 

where  nQj  and  nj j  (£,t)  are  the  equilibrium  and  perturbed  densities,  respective¬ 
ly.  In  (12)  v>m,  where  j=m  denotes  the  Maxwellian  electrons,  is  approximately 
the  elastic  collision  frequency,  and  vh  is  the  sum  of  the  elastic  and  vibra¬ 
tional  excitation  collision  frequency  for  the  suprathermal  electrons.  The 
important  point  to  note  is  that  the  Maxwellian,  not  the  suprathermal,  electrons 
are  effective  in  damping  the  waves,  and  that  (12)  is  a  good  approximate  operator 
for  the  study  of  damping  effects. 

Assuming  perturbations  of  the  form  exp( il<.£-iu-t)  we  solve  Eq.  (10)  for  fjj 
and  then,  with  the  aid  of  Eq .  (11),  obtain  the  dispersion  relation 
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Here  fl=qB0/mc  is  the  electron  cyclotron  frequency,  u>p  *  4nq2  n0/m,  where  nG  = 

^  noj  *  So j  =  foj/noj  »  P j  =  Ooj/no  ,  k2  *  k2  +  k2  »  and  other  quantities  have 

their  usual  meanings.  For  the  Maxwellian  electrons  with  density  n^  and  tem¬ 
perature  Te,  the  velocity  space  Integrations  can  be  carried  out  and  Eq.  (13) 
becomes 

•4m  k«  ,  .  ..  “  2ann2<J2mki(aH-ivm) 

“  u.<a.+ivm)k2  Lo*XP{~b)  "  n-1  a>[Ca.+ivm>2-T12n2]k2 

+  (<U2/«2)  (Nh/Dj,)  -  0  .  (16 

Here  an  =  (In/b)  exp(-b) , In(b)  is  the  modified  Bessel  function  of  the  first 
kind,  b  >=  k2Te/mft2  ,  *  4irq2  nom/m.  We  have  assumed  k2<.<.k^  ,  and  (u)-nft)/ 

[  |k|  |(2Te/  m  )1/2]»1  for  all  n,  so  that  Landau  damping  and  cyclotron  damping  can 
be  neglected. 


1.1. 2.1  Upper  Hybrid  Instability 

In  the  lower  ionosphere  we  may  assume  b<l,  and  ft2/u>241.  Without  the  supra 

thermal  electron  term,  one  of  the  roots  (u>-u)r+iuic  ,  u>c«wr)  of  the  dispersion 
equation  (16)  is 


(17) 


wr  =  (a.2  +  n2)1^2 

r  v  pm 

which  is  the  upper  hybrid  frequency,  and 

(0C  =  -  [(<^n  +  02)/2«j2]  Vm  ,  (18) 


is  the  collisional  damping  rate. 


The  suprathermal  electrons  interact  resonantly  with  the  upper  hybrid  wave 
and  drive  it  unstable,  the  resonance  condition  being  u.r-k||V||=  n ft.  Since  ^ 
nom,  the  growth  rate  Y  can  be  assumed  to  be  small  compared  to  o.r  and  is  given  by 
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The  collisions  of  the  suprathermal  electrons  will  somewhat  detune  their  reso¬ 
nance  with  the  wave  and  reduce  the  growth  rate.  This  damping  rate  can  be  shown 
to  be  “(noh/non,)  vj,.  So,  the  significant  damping  is  that  due  to  the  Maxwellian 
electrons  and  is  given  by  Eq .  (18).  The  explicit  expression  for  Y  in  the  energy 
representation  using  G0j,(E)  =  (2^2  E^/2/m^/2)  8oh(v)  *  80  that  4x/dE  G0h  =  1,  is 
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Here  E0=mvo/2  is  the  peak  energy  (-5  eV)  of  the  suprathermal  electrons,  Kj|  = 

k||v0/ft  »  En  “  Eo[((‘-r/n_n)/Kll  l2  *  and  an  =  (kiv0/n)  [  (E-En)/E0] 1^2 .  We  notice 
that  Y  depends  on  the  parameters:  noh/nom>  Wp/0,  kj ,  and  kx ,  which  are 
altitude  dependent.  At  a  given  altitude  the  unstable  spectrum  is  given  by  the 

condition  J2  .  -  J2  ,  ^0. 

n+1  n-1 

Numerical  values  of  Y  at  various  altitudes  have  been  obtained  using  the 
energy  distributions  of  Jasperse*,  and  the  results  are  given  in  Table  1.  The 
upper  hybrid  instability  can  be  excited  at  an  altitude  as  low  as  125  km,  below 
which  it  is  stabilized  by  the  collisions.  Above  170  km,  the  damping  effects  of 
electron-electron  and  electron-ion  collisions  must  be  considered. 


■  —  *  -  t  «  *  _  *  < 
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TABLE  1.  Values  of  Upper  Hybrid  Frequency  (o^) ,  Suprathermal  Electron 
Abundance  (a).  Elastic  Collision  Frequency  of  Maxwellian 


I. 1.2. 2  Electron  Cyclotron  Instability 


Here  we  consider  waves  with  kj=0,  and  the  instability  that  is  found  is  of 
the  nonresonant  type  with  u)r=nft.  In  this  case  we  can  write 
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(22) 


where  o  =  noh/nom  =  8h,  and  v2  =  v2  +  v2 .  We  first  solve  Eq.  (16)  with  kH=0  in 
the  limit  vm=0.  We  substitute  a)  =  nft(l+X).  where  n>2  and  ,  and  obtain 
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where  Ki  =  kiv0/ft,  o)2  K  cl)2  ,  and 

p  pm 

I  =  Att Eq  Jj  Goh(E)  J2n  [2K1(E/E0)1/2]  (2A) 

in  the  energy  representation.  We  have  neglected  terms  of  the  order  of  (v^/fi)  o 
or  smaller.  It  is  readily  seen  that  x  will  have  a  positive  imaginary  part 
(meaning  instability)  if  1*0,  |l | >  an  K2/o,  and 
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Condition  (25)  determines  the  range  of  values  of  oJ2/ft2  for  which  the  n-th 

P 

harmonic  can  be  unstable.  The  maximum  values  of  Imx  is  obtained  when  Wp/&2  = 
(n2-l)2ai«  In  this  case  x  Is  given  by 
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The  second  term  within  the  parentheses  represents  the  finite-b  stabilizing 
effect.  In  the  presence  of  collisions  (vm*0),  the  complete  solution  of  the 
dispersion  relation  is 


a) 
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(n2+l)vm 
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(27) 


where  Y  =  nft(Imx)  and  the  additional  term  is  the  collisional  damping  rate. 

The  electron  cyclotron  instability  has  a  larger  growth  rate  than  that  of 
the  upper  hybrid,  and  an  unstable  harmonic  is  localized  within  a  narow  range 
(~2  km)  around  the  altitude  at  which  It  has  the  maximum  growth  rate.  Results  of 
numerical  calculations  are  summarised  in  Table  2.  At  altitudes  100  km  and  below, 
the  mode  Is  stabilized  by  the  electron-neutral  collisions  while  above  130  km  it 
is  stabilized  by  the  finite-b  effect. 


1.2.2.  Quasillnear  Evolution  of  the  Waves  and  of  the  Photoelectron  Distribution 


We  have  studied  the  upper  hybrid  waves  analytically  in  somewhat  details. 
Excitation  of  these  waves  requires  resonant  interaction  between  the  suprathermal 
electrons  and  the  wave,  and  the  resonant  diffusion,  in  velocity  space,  of  the 
equilibrium  suprathermal  electrons  in  response  to  the  unstable  waves  is  governed 
by  the  following  equation: 
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where  the  time  evolution  of  the  field  energy  density  <f(k,t)  is  given  by 

57  (f(k,t)  =  2  Im  oj(k,t)  £(k,t)  ,  (29) 

and  and  are  given  by  Eqs.  (14)  and  (15)  with  j=h,  respectively. 

Equation  (28)  may  be  derived  by  following  the  standard  method  of  the  quasilinear 
theory.  The  superscript  R  on  g0h  is  to  remind  us  that  the  equation  is  valid  in 
the  resonant  region. 

An  H-theorem  may  be  demonstrated  from  Eq.  (28)  by  multiplying  by  g^  and 
integrating  over  v^.  This  gives 

1  d  fa  r2  .  .  8ir2e2  r  r  £(k,t)  , 

2  3Fjdi  eoh  ^  •  -  —jr  „„1  JdiJdt  —2 


« d2„  <^> 


r  r  (s£L  _i_ 


+  k|  a^>  eot,)2 


'  £  k2  [I"  «■'*>!*  5^  k*  377>  4  •  <3°> 


Referring  to  Eqs.  (14)  and  (15)  we  see  that 
t  /.X  x  _  ^  noh  f. 


,  u*-  uoh  ^  r  o  kivi 

Tm(Nh/Dh)  -  -it  J  dv  ( — — )  5  (u)r-k||V|,-nS) 


,nQ  3  3  . 

x  (n  5^1  k|1  37^  8oh 


Comparing  the  two  terms  inside  the  curly  brackets  on  the  right-hand  side  of 
Eq.  (30)  we  conclude  that  the  right-hand  side  of  Eq .  (30)  is  negative.  Equation 
(30)  then  states  that  the  time  rate  of  change  of  a  positive  quantity  is  non¬ 
positive.  Consequently,  the  time-asymptotic  state  must  be  the  one  for  which  the 
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integrand  of  the  velocity  integration  on  the  right-hand  side  of  Eq .  (30) 
vanishes.  That  is, 


/; 


dki  <S(kifki  = 
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c 


where  g^  (vy®)  =  gkh  (.)!.» t=co)>  an&  Nh(“)»  Dh(°°)  are  to  be  evaluated  at  k u  = 

(u)r-nft)/vn •  It  is  expected  that  the  field  fluctuations  will  grow  in  time  and 
then  saturate  at  a  finite  amplitude  level  in  the  time-asymptotic  state.  In 


other  words,  £(t=®°)*0.  Then,  according  to  Eq.  (32),  g8^(v,°°)  Is  given  by 


[<£»  J_  +  ^  J_)  gR  (v.,l2 

+  3V|, '  8oh 


V1  3v^  v  ||  3  v  ||  6oh 
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At  this  point  we  postulate  that  which  is  initially  isotropic  in  velocity 


space  remains  so  at  all  times.  Then 
nft  3  (dr-nfl  3 


(—  3 - *" 
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where  v  =  v 


,  and  (33)  becomes 


±  gR  (v,-)]2  .  JS  iH  %  .  —  gR 
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Referring  to  Eq .  (19),  and  using  (34)  we  find 


2«2 


Im  ;  Nh(“>)/Dh(°°)  ]  = - 7—  wr  y(“) 


c 


where  Y(°°)  is  the  time-asymptotic  collisionless  growth  rate.  Next,  we  assume 
that  Y  (°°)  =  |a)c  |  so  that  the  fluctuations  do  not  grow  any  more.  Here  uc  is  the 
collisional  damping  rate  and  may  be  taken  as  given  by  Eq.  (18).  Using  this 
assumption  in  Eq.  (35)  we  find 
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v  g 


R  —  g 
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(37) 


which  can  be  solved  to  yield 


8oh^V,°°)  =  C  exP(~av2/2)  > 
in  the  resonant  region  of  the  velocity  space,  where 


no  vmvh  ,  2 

a  =  -  — —  kf  , 

noh  cn4  i 

P 

and  where  the  constant  C  has  to  be  determined  by  matching  this  solution  to  the 
rest  of  the  distribution  function.  Changing  over  to  the  energy  representation 
we  find 

Goh(E’“)  *  £l/2  exP(~a'E/E°)  »  (38) 


where  E0  is  the  peak  energy  of  the  supra thermal  electrons  and 


>» 

a 


1  J^O  a4 

2  noh  fl2  Wp 


(39) 


with  Ki  s  kj.v0/fi.  Using  the  values  of  the  parameters  involved,  a'  is  estimated 
to  be  ~10_2.  The  expression  (38)  is  valid  for  E>En,  where  En  has  been  defined 
before,  corresponding  to  the  resonant  region  of  the  velocity  space.  We  infer 
that  in  the  resonant  region  the  suprathermal  electron  distribution  function, 
which  initially  had  a  positive  slope,  will  have  a  small  negative  slope  in  the 
time-asymptotic  state  due  to  the  wave-induced  diffusion  process.  Of  course, 
when  the  photoelectron  source  function,  which  is  responsible  for  the  bump  in  the 
equilibrium  distribution  function  at  t=0 ,  is  retained  in  the  analysis  the 
distribution  function  in  the  time-asymptotic  state  may  have  a  small  bump  in  the 
2-6  eV  region.  For  this  a  numerical  analysis  is  required.  Our  conclusion  is 
that  the  experimentally  measured  photoelectron  distribution  in  the  2-6  eV  energy 
range  can  be  well  explained  if  waves  and  wave-induced  diffusion  process  are 
included  in  the  theory  of  Jasperse^-. 


1-17 


II.  MODELING  OF  DAYTIME  AIRGLOW  AND  ELECTRON  DENSITY  PROFILE  IN  THE 

MIDLATITUDE  IONOSPHERE 

Because  of  its  influence  on  many  Air  Force  communication  and  surveillance 
systems,  it  is  of  great  importance  to  know  the  electron  density  profile  (EDP)  of 
the  ionosphere.  In  circumstances  where  ionospheric  sounders  and  other  direct 
means  cannot  be  used  to  determine  this  information,  some  means  of  modeling  the 
EDP  or  determining  it  by  remote  sensing  must  be  found. 

One  useful  diagnostic  of  conditions  within  a  plasma  is  the  intensity  of 
optical  emission  features  at  different  wavelengths.  In  the  daytime  ionosphere, 
electrons  released  by  photoionization  excite  atoms  and  molecules,  which  then 
emit  photons  at  characteristic  frequencies  when  they  return  to  the  ground  state. 
Locally  the  emission  rate  depends  on  the  density  of  the  emitting  species  as  well 
as  the  spectrum  of  the  photoelectron  flux.  The  overall  emission  measured  along 
a  line  of  sight  will  depend  upon  the  profiles  of  these  quantities. 

In  an  effort  to  evaluate  the  practicality  of  using  optical  emissions  as  a 
remote  sensor  of  the  EDP,  we  have  developed  the  capability  of  modeling  the 
airglow  emission  of  several  features  in  the  ultraviolet  -  the  LBH  bands  and  the 
spectral  lines  at  1356A  -  emitted  by  the  daytime,  midlatitude  ionosphere.  For 
the  daytime  no  method  is  yet  known  for  determining  the  EDP  directly  from  optical 
emissions;  instead  an  Indirect  method,  using  optical  emissions  to  constrain  the 
parameters  of  a  model  of  the  EDP,  will  likely  be  necessary.  In  this  report  we 
describe  a  series  of  case  studies,  modeling  the  EDP  using  first-principle  calcu¬ 
lations  and  comparing  the  results  to  a  variety  of  direct  measurements.  This 
series  is  aimed  at  evaluating  how  well  such  ab  initio  calculations  can  model  the 
EDP  and  at  determining  which  parameters  of  the  model  most  sensitively  control 
the  EDP.  For  cases  in  which  simultaneous  optical  emission  measurements  are 
available,  we  evaluate  the  emissions  predicted  by  our  model  EDP  and  compare  with 
the  observed  emissions. 

In  our  ab  initio  calculations,  the  electron  and  ion  densities  are  found  by 
solving  the  continuity  equations  for  all  species  simultaneously.  In  the  bottom- 
side  ionosphere  transport  can  be  ignored,  and  in  the  steady-state  approximation 
the  continuity  equation  for  a  species  simplifies  to  a  balance  between  its  pro¬ 
duction  and  loss  rates.  The  four  most  important  ion  species  -  0+,  N2+,  02+,  and 
N0+  -  are  included  in  our  model. 
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Above  the  EDP  peak,  the  full  time-dependent  continuity  equation  including 
transport  must  be  solved,  although  a  simplification  is  achieved  because  only  the 
dominant  ion  species,  0+,  need  be  followed.  We  have  modified  an  existing 
transport  code  so  that  it  can  run  concurrently  with  our  bottomside  code.  The 
result  is  a  calculated  EDP  from  100  to  1000  km. 

To  calculate  the  intensity  of  VUV  and  UV  daytime  emissions,  we  begin  with 
the  photoelectron  flux,  F(E,z),  calculated  as  a  function  of  energy  and  altitude. 
The  volume  emission  rate  at  altitude  z  is  given  by 

Si(z)  =  n(z)  J  F(E,z)  o(E)  dE 

where  n(z)  is  the  density  of  the  emitting  species  and  o  is  the  corresponding 
cross  section.  Integration  of  S^(z)  along  the  line  of  sight  taking  into  account 
absorption  gives  us  the  column  emmission  rate. 

The  key  parameters  of  these  models  are:  1)  the  solar  EUV  flux  2)  the 
neutral  atmosphere  (N2,  O2  and  0  densities),  3)  the  neutral  wind  and  electric 
fields,  and  4)  the  temperatures  (Te,  T^  and  Tn)-  Although  the  latter  two  are 
strictly  parameters  of  the  transport  model  for  the  EDP  above  the  peak,  the  first 
two  play  a  major  role  in  determining  the  EDP  throughout  the  daytime  ionosphere. 

It  is  unlikely  that  all  the  parameters  needed  to  calculate  an  ab  initio 
model  of  the  EDP  will  be  known  in  any  particular  case.  To  judge  how  much  we 
could  rely  on  estimates  of  unknown  quantities,  and  how  well  our  codes  would  per¬ 
form,  we  conducted  a  series  of  tests  under  such  circumstances,  in  which  we  could 
compare  the  results  of  our  ab  initio  calculations  directly  with  measured 
electron  density  profiles. 

(1)  Case  One:  White  Sands  Missile  Range,  August  23,  1972,  Heroux 

et  al . ,  (1974). 

This  case  involves  rocket  data  for  the  solar  EUV  flux  coincident  with 
ground  based  ionosonde  data.  This  represents  a  partial  control  on  the  produc¬ 
tion  rate  and  should  lead  to  good  bottomside  EDP  agreement.  In  Figure  II-l,  we 
show  a  plot  of  the  bottomside  electron  density  profile  from  the  ionosonde 
measurement  and  a  graph  of  our  hybrid  calculation  for  the  electron  density  pro¬ 
file.  It  should  be  noted  that  the  estimated  errors  for  the  ionosonde  measure¬ 
ment  were  130  km  near  120  km  and  15  km  elsewhere.  For  this  case  a  Jacchia 
(1977)  neutral  atmosphere  with  the  appropriate  parameters  for  the  given  day  was 
used  in  the  calculation  in  addition  to  the  measured  solar  EUV  flux.  Typical 


daytime  models  were  used  for  the  electron  temperature  and  neutral  winds.  With 
observational  control  on  the  solar  flux  one  would  expect  good  EDP  agreement  in 
and  below  the  layer.  The  agreement  in  the  Fj  layer  is  in  fact  within  15%, 
while  the  agreement  in  the  E  region  is  within  the  error  bars  on  the  ionosonde 
EDP.  Unfortunately  the  unusually  large  error  bar  on  the  ionosonde  measurement 
near  120  km  renders  the  E  region  agreement  of  limited  quantitative  significance. 
The  25  km  or  30%  difference  in  the  F2  peak  of  the  EDP  is  consistent  with  any 
other  state-of-the-art  ab  initio  model  calculation,  and  merely  reinforces  our 
intent  to  pursue  the  study  of  the  transport  term  which  is  in  all  likelihood  the 
source  of  this  difference. 

(2)  Case  Two:  White  Sands  Missile  Range,  August  14,  1979 

Rocket  data  for  the  solar  EUV  flux  coincident  with  ground-based  ionosonde 
data  were  again  available  for  this  case.  Thus,  we  could  use  the  measu .d  solar 
flux,  along  with  the  appropriate  Jacchia  atmosphere  model  and  the  neutral  wind 
and  temperature  models  used  in  the  previous  case  to  calculate  the  ab  initio 
model.  The  results,  shown  in  Figure  II-2,  are  similar  to  those  of  the  previous 
case,  although  the  accuracy  of  our  calculated  peak  altitude  happens  to  be  better 
in  this  case,  with  the  altitude  of  the  calculated  peak  falling  within  5  km  of 
the  measured  position. 

(3)  Case  Three:  ISIS  2  topside  sounder.  May  23,  1972 

The  data  considered  here  come  from  an  ISIS  2  pass  at  8:28  UT  at  60°  north 
latitude  and  38°  east  longitude.  The  data  were  taken  from  iso-density  contours 
in  a  sample  collection  of  ISIS-2  observations  (Klumpar,  1980).  Due  to  reported 
difficulties  in  interpretation  of  sounder  data,  we  have  shown  the  data  with  150 
km  error  bars.  Figure  II-3  shows  two  theoretical  curves  with  the  topside  sounder 
data.  The  solid  curve  shows  the  calculated  profile  using  our  best  a  priori 
estimates  for  the  input  quantities.  We  can  see  that  the  ab  initio  calcula¬ 
ted  profile  is  60%  lower  than  the  data  but  agrees  in  shape.  The  dotted  curve  is 
a  scaled  version  of  the  calculated  profile  fitted  to  the  data.  Had  this  been  a 
nighttime  case,  we  would  expect  direct  VUV  emissions  as  discussed  in  the  next 
section  to  provide  a  good  estimate  of  the  peak  of  the  EDP,  to  which  the  calcu¬ 
lated  topside  EDP  shape  would  then  be  scaled.  This  would  lead  to  quite  satis¬ 
factory  agreement  for  this  case.  For  the  daytime  case,  the  dashed  curve  would 
be  essentially  duplicated  above  350  km  by  use  of  a  stronger  equatorward  wind  in 
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the  model  (as  might  follow  at  these  latitudes  from  a  high  latitude  heating 
event).  The  excellent  agreement  of  the  shape  of  the  model  and  observed  EDP 
satisfies  necessary  conditions  for  the  merit  of  this  EDP  scheme.  However,  some 
additional  observable  must  be  added  to  scale  the  ab  initio  calculated  EDP  before 
we  can  look  for  the  good  agreement  shown  by  the  dashed  curve. 

(4)  Case  Four:  Millstone  Hill,  incoherent  scatter  radar,  April  8,  1978 

This  is  another  case  in  which  the  solar  EUV  flux  and  other  input  parameters 
for  our  model  are  unknown  and  we  must  attempt  to  estimate  them.  Using  typical 
values  for  these  parameters  we  calculated  the  EDP.  Figure  II-4  shows  a  com¬ 
parison  with  the  EDP  derived  from  an  incoherent  scatter  measurement.  In  this 
case  the  predicted  peak  is  around  50  km  above  the  measured  peak,  and  densities 
can  be  as  much  as  a  factor  of  three  in  error. 

Our  conclusion  from  these  four  cases  is  that  when  we  can  specify  the  input 
data  necessary  for  our  EDP  model,  the  model  does  well  in  fitting  the  actual  EDP. 
When  all  the  necessary  data  is  not  available,  our  model  does  a  poorer  job. 

We  now  turn  to  cases  in  which  optical  emission  data  are  available. 

(5)  Case  Five:  S3-4 ,  Rev  373,  April  8,  1978 

In  this  case  we  had  detailed  optical  emission  data  from  a  nadir-looking  VUV 
spectrometer  flying  onboard  DoD  satellite  S3-4,^  but  no  simultaneous  solar  EUV 
flux  or  EDP  measurements.  To  calculate  airglow  emission  intensities  we  need,  in 
addition  to  the  parameters  of  our  EDP  model,  the  cross  sections  for  the 
radiation  processes.  Because  of  the  difficulty  of  laboratory  mesurements,  some 
cross  sections  are  still  not  well  known.  In  Figure  II-5  we  show  airglow  calcula¬ 
tions  from  our  ab  initio  model  for  features  at  nine  wavelengths,  along  with  com¬ 
parisons  with  the  S3-4  data.  We  see  that  agreement  is  good  at  low  wavelengths, 
around  the  feature  at  1356A,  but  that  the  agreement  noticeably  worsens  at  longer 
wavelengths  where  scattered  light  may  be  contaminating  the  data.  The  large 
error  bars  on  the  data  points  are  due  to  the  statistical  uncertainty  resulting 
from  the  low  count  rates  in  the  data. 

In  Figure  I 1—6  we  track  the  optical  emissions  at  1356A  as  the  satellite 
moves  along  its  orbit  in  the  noon-midnight  plane,  passing  through  midday  in  the 
northern  hemisphere.  We  find  that  our  completely  ab  initio  calculation  does  a 
remarkable  job  in  reproducing  the  observed  emissions  produced  under  a  wide  range 
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of  conditions.  A  similar  plot  for  the  feature  at  1383A  (Figure  II-7),  shows  our 
theory  consistently  low  compared  to  the  data. 

(6)  Case  Six:  Hilat  Satellite  July  1983 

In  this  case  we  attempted  to  model  the  1356A  airglow  observed  from  the 
Hilat  satellite.  On  the  date  of  the  observations  the  satellite  passed  overhead 
near  a  ground  based  ionosonde  at  Millstone  Hill  making  near-coincident  measure¬ 
ments  of  the  EDP  available. 

The  optical-emission  measurements  were  made  by  the  AIM  sensor  onboard  Hilat 
at  22  hr  15  min  10  sec  UT  on  Rev.  219.  At  that  time,  while  crossing  the  lati¬ 
tude  of  Millstone  Hill,  the  satellite  passed  16.4°  east  of  Millstone  Hill.  The 
AIM  sensor  scans  through  135°  along  a  path  perpendicular  to  the  satellite  orbit, 
in  this  case  beginning  its  sweep  looking  east  into  darkness  and  ending  looking 
west  into  the  bright  limb,  Figure  II-8.  There  is  some  uncertainty  in  the  orien¬ 
tation  of  the  satellite,  because  in  this  rev.  it  was  still  oscillating  on  all 
three  axes  with  the  amplitude  of  the  oscillation  unknown.  If  we  assume  that  the 
sensor  was  pointed  in  the  nadir  direction  at  the  midpoint  of  the  sweep,  then  the 
observed  direction  of  the  bright  limb  is  displaced  by  one  data  bin,  5.64°,  from 
its  expected  location.  We  assume  that  this  gives  us  an  estimate  of  the 
satellite  roll,  and  we  assign  error  bars  of  this  magnitude  to  the  directions  of 
the  calculated  optical  emissions.  The  look  direction  from  the  satellite  to 
Millstone  Hill  is  53.62°W,  placing  Millstone  Hill  in  the  foot  of  the  bright 
limb.  An  uncertainty  of  5.6°  in  the  look  direction  towards  Millstone  forces  an 
uncertainty  in  the  1356A  flux  of  the  order  of  25%. 

Our  modeling  began  with  a  nadir  calculation  at  the  satellite  location,  with 
a  solar  zenith  angle  (SZA)  of  81°.  At  this  high  SZA,  the  plane  parallel 
approximation  for  the  solar  flux  transport  breaks  down,  resulting  in  an 
underestimation  of  the  airglow  emission  rates.  To  get  an  upper  bound  for  the 
airglow,  a  second  nadir  calculation  was  performed  using  an  artificial  value  of 
the  SZA,  the  angle  whose  secant  approximates  the  Chapman  function  at  81°,  that 
is,  79°.  The  results  of  the  two  calculations  are  plotted  in  Figure  II-9  at 
pixel  160  (different  pixels,  are  different  look  directions)  we  find  that  they 
are,  respectively,  1.2  and  1.5  times  the  observed  1356A  airglow  intensity  In 
that  direction. 
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To  calculate  the  column  emission  rate  at  some  angle  to  nadir  a  plane 
parallel  approximation  for  the  atmosphere  was  used.  Calculations  at  10,  20,  30, 
40  and  53.62°  (the  look  direction  to  Millstone  Hill)  are  presented  in  Figure 
II-9.  We  see  that  the  calculated  emission  profile  reflects  the  profile  of  the 
data.  At  53.62°,  however,  the  change  in  SZA  along  the  look  direction  was  large 
enough  to  call  the  plane-parallel  approximation  into  question.  At  that  point  we 
performed  two  airglow  calculations  to  give  bounds  to  the  emission.  One  was 
based  on  the  volume  emission  rate  below  the  satellite,  which  with  its  higher  SZA 
gives  a  lower  bound  to  the  airglow.  The  other  calculation,  based  on  the  volume 
emission  rate  above  Millstone,  gives  an  upper  bound.  Figure  II-C  shows  the  two 
results  differing  by  a  factor  of  two. 

We  were  Interested  in  this  particular  sweep  because  a  ground-based 
ionosonde  had  measured  the  EDP  just  a  half  hour  before  the  satellite  fly-over. 
While  not  an  ideal  coincidence  either  spatially  or  temporally,  this  did  give  us 
an  opportunity  to  compare  our  ab-initio  calculations  for  the  EDP  and  the  airglow 
with  observations  of  both  quantities. 

Comparison  of  the  measured  EDP  with  the  ab  initio  EDP  calculated  for  the 
airglow  theory  presented  above  shows  that  the  ab  initio  EDP  is  consistently  low 
(Figure  TI-10).  We  are  currently  working  to  find  a  method  of  adjusting  the 
input  parameters  of  our  EDP  calculation  which  can  bring  the  calculated  EDP  into 
line  with  the  measured  EDP  without  destroying  the  fit  of  the  calculated  airglow 
eml-sion  to  the  observed  airglow.  One  such  adjusted  model,  shown  in  Figure 
1 1  —  1 0 ,  was  made  by  changing  the  neutral  atmosphere  to  fit  the  model  EDP  to  the 
measurements  at  the  F2  peak.  Such  an  adjustment,  however,  will  also  increase 
the  already  high  optical  emission  rate. 

The  conclusion  of  this  study  then  is  that  our  ab  initio  calculations  can 
perform  fairly  well  when  crucial  input  data  can  be  accurately  specified. 

Airglow  measurements  may  be  helpful  in  constraining  parameters  when  the  input 
data  is  inadequate.  The  future  direction  of  this  work  should  be  aimed  at 
establishing  the  correlations  between  the  airglow  and  the  EDP  parameters  so  that 
the  airglow  measurements  can  be  used  for  this  purpose. 
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III.  PARTICLE  PRECIPITATION:  MAGNETOSPHERE-IONOSPHERE  COUPLING 


The  nighttime  ionosphere  at  high  latitudes,  i.e.,  in  the  auroral  zone,  Is 
created  by  energetic  particles  (both  electrons  and  ions)  precipitating  from  the 
magnetosphere.  Calculations  of  energy  deposition  and  ionization  in  the 
atmosphere^- *2  require  that  the  incident  flux  of  particles  precipitating  from 
the  magnetosphere  be  specified.  We  have  addressed  this  question  with  a  study  of 
the  pitch  angle  diffusion  process  in  the  magnetosphere  which  causes  particles  to 
precipitate.  The  results  of  this  study  are  summarized  in  subsection  111.1,  below 
We  were  able  to  calculate  the  loss-cone  population  in  the  magnetosphere,  as  a 
function  of  position,  when  the  processes  of  precipitation  and  diffusion  compete 
against  each  other;  the  results  are  illustrated  using  data  from  the  diffuse 
aurora. 

In  addition  to  ionization,  the  precipitating  particles  can  also  trigger 
plasma  instabilities  in  the  ionosphere,  which  can  cause  anomalous  heating,  irregu 
larities,  and  other  phenomena  which  interfere  with  communications.  One  of  the 
sure  indications  of  such  activity  in  the  aurora  zone  is  the  detection  of  ion 
conics:  fluxes  of  energetic  ions  streaming  out  of  the  ionosphere  along  auroral 

field  lines.  We  have  studied  both  the  means  by  which  turbulence  can  be  created 
by  precipitating  electrons  and  the  way  in  which  the  turbulence  can  accelerate 
ionospheric  ions.  This  work  is  summarized  in  subsection  III. 2,  below.  In  an 
empirical  study,  using  a  Monte  Carlo  model  for  wave  particle  interaction,  we 
were  able  to  show  that  the  commonly  observed  intense  lower-hybrid  turbulence  in 
the  suprauroral  region  is  capable  of  accelerating  ions  to  the  energies  commonly 
observed  in  ion  conics.  Then,  using  a  plasma  simulation,  we  demonstrated  that  a 
precipitating  electron  flux  modeled  after  the  flux  observed  in  a  discrete 
auroral  arc  could  excite  lower-hybrid  turbulence  and  cause  significant  ion  acce- 


Ill .1 .  The  Loss-Cone  Population  In  the  Magnetosphere 

Pitch-angle  scattering  has  long  been  recognized  as  one  of  the  important 

3 

consequences  of  wave-particle  interaction  in  the  magnetosphere  .  Its  action 
continuously  replenishes  the  flux  of  particles  in  the  loss  cone  in  velocity 
space,  which  gain  access  to  the  dense  atmosphere  and  precipitate  out  of  the 
magnetosphere^.  Indeed,  for  the  continuous  (diffuse)  aurora  a  common  approxi¬ 
mation  is  to  assume  that  scattering  is  strong  enough  to  maintain  an  isotropic 
particle  distribution  in  the  face  of  such  anisotropic  particle  loss^*.  Detailed 
satellite  observations^,  however,  show  that  electrons  at  keV  energies  and  above 
are  often  anisotropic.  At  the  low  altitude  of  the  ISIS  satellite,  the  upward 
flowing  loss  cone,  containing  the  few  electrons  backscattered  from  the 
atmosphere,  is  always  strongly  depleted.  It  is  more  interesting  to  find  that 
the  downward  flowing  loss  cone,  containing  precipitating  electrons,  often  shows 
a  moderate  depletion,  too.  While  pitch-angle  scattering  has  attempted  to  fill 
in  the  strongly  depleted  loss  cone  created  at  the  other  footpoint,  its  job  is 
still  unfinished  when  the  particles  reach  the  end  of  the  field  line  and  precipi¬ 
tate  . 


The  effect  of  pitch-angle  scattering  by  wave-particle  interaction  can  be 
described  by  the  quasilinear  diffusion  equation 
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In  this  equation,  f  is  the  distribution  function,  E  and  B  are  the  static  fields 
and  D  is  the  quasilinear  velocity  diffusion  coefficient  due  to  fluctuating 
fields.  (In  a  self-consistent  treatment  E  would  be  determined  from  the  particle 
populations  via  Poisson’s  equation  or  the  condition  of  quasineutrality^  and  the 
diffusion  coefficient  would  be  determined  from  the  turbulence  arising  from  the 
instability  of  features  in  the  particle  distributions^.  For  the  sake  of  simpli¬ 
city,  we  regard  these  fields  as  given,  imposed  on  the  particles.) 

When  the  spatial  variation  of  the  distribution  function  is  weak,  a  fruitful 
simplification  of  this  equation  can  be  made  by  averaging  over  the  bounce  motion 

of  a  particle  .  The  resulting  bounce-averaged  diffusion  equation  has  seen  wide 
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application  in  the  magnetosphere. 

Far  outside  the  loss  cone  the  bounce-averaged  diffusion  equation  can  be 
justified  because  the  spatial  gradients  are  small.  Inside  the  loss  cone, 
however,  the  bounce-averaging  approximations  fail.  When  particles  are  removed 
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from  the  loss  cone,  sharp  gradients  in  the  distribution  function  are  created  at 
its  edge.  These  gradients  magnify  the  diffusion  term  so  that  significant  change 
in  the  distribution  function  can  occur  over  a  time  short  compared  to  a  bounce 
period.  This  implies  significant  variation  of  the  distribution  function  along 
the  orbit,  which  the  bounce-averaged  equation  cannot  describe. 

To  calculate  the  spatially  varying  distribution  function,  we  must  return  to 
the  full  equation  (1).  Numerical  solutions  of  (1)  for  the  distribution  function 
over  all  velocity  space  have  been  calculated  by  Davidson‘S •  But  the  narrowness 
of  the  loss  cone  may  be  turned  to  our  advantage,  for  it  implies  that  we  need  to 
solve  the  spatially  dependent  problem  only  for  that  small  region  of  velocity 
space.  A  boundary-layer  treatment  along  this  line  has  been  developed  for  the 
loss  cone  in  a  mirror  machine  plasma^2  and  for  the  stellar  distribution  around  a 
black  holeS .  Basically  one  solves  Eq .  (1)  in  a  simplified  form  for  the  loss 
cone,  taking  advantage  of  the  smallness  of  the  velocity-space  region  over  which 
the  solution  is  to  apply,  and  matches  that  solution  to  a  solution  of  the  bounce 
averaged  equation  which  applies  outside  the  loss  cone.  We  will  use  this  method 
to  describe  the  spatially  dependent  distribution  function  in  the  loss  cone  for 
the  geomagnetic  field.  We  will  then  apply  the  results  to  interpret  particle 
data  from  a  satellite  pass  through  the  continuous  aurora.  Readers  interested  in 
more  detail  than  is  presented  here  are  referred  to  Retterer  et  al.^. 

111.1.1  The  Boundary -Layer  Equation 

The  kinetic  equation  (1)  describes  both  the  random  changes  in  velocity  due 
to  scattering  and  the  systematic  changes  with  position  due  to  the  mean  fields. 

A  simplification  of  the  equation  results  when  a  change  is  made  to  variables 
incorporating  the  systematic  changes,  viz .  particle  energy  and  magnetic  moment, 
or  the  magnitude  of  the  velocity,  vQ,  and  the  pitch  angle,  a0,  of  the  particle 
when  it  crosses  the  equatorial  plane.  Assuming  that  all  other  velocity  gra¬ 
dients  will  be  smaller,  we  retain  only  the  pitch-angle  gradients  in  the  dif¬ 
fusion  term.  We  will  look  for  a  quasisteady  state  solution  for  f:  after 
possible  initial  transients  the  only  time  variation  we  expect  is  the  slow  decay 
of  the  population  as  the  result  of  precipitation  through  the  loss  cone.  With 
these  approximations,  the  kinetic  equation  becomes 
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where  J  is  the  Jacobian  for  the  transformation  from  local  velocity  v  to 
equatorial  velocity  v0,  and  D0  is  the  familiar  equatorial-pitch-angle  diffusion 
coefficient  (with  units  of  sec-1).  We  simplify  the  equation  further  by  intro¬ 
ducing  a  scattering  depth  t,  normalized  to  run  from  0  to  1  as  s  runs  from  one 
footpoint  of  the  field  line  -s^,  to  the  other  at  +sm: 
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Tb  Is  the  bounce  period,  and  <D0>b  is  the  bounce  averaged  diffusion  coef¬ 
ficient.  Introducing  the  pitch-angle  variable  u 

u  =  sin2a0/2X  (5) 

and  using  the  small-angle  approxima t Ion  slna0*aQ,  we  reduce  the  kinetic  equation 
to  Its  final  form: 
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1 1 1 . 1 . 2  The  Pi  t  chrAn_gl  e_  Distribution 

We  solved  Eq .  (b)  with  the  appropriate  boundary  conditions  to  find  the 
distribution  function  as  a  function  of  pitch  angle  (u)  and  position  (t).  The 
solution  lias  two  parameters.  The  first  is  u£,  the  value  of  u  at  the  edge  of  the 
loss  cone 

U£  =  s  1  n2ot0£  /  2  Tb<D0>b  .  (7) 

This  determines  the  degree  of  anisotropy  within  the  loss  cone;  the  larger  u£  is, 
the  larger  is  the  size  of  the  loss  cone  compared  to  the  mean  angle  through  which 
a  particle  is  scattered  during  a  bounce  period  and  thus  the  greater  is  the 
depletion  within  the  loss  cone.  The  second  parameter  is  an  arbitrary  nor¬ 
malization  factor,  which  Is  determined  (when  needed)  by  matching  the  boundary- 
layer  solution  to  the  bounce-averaged  solution  outside  the  loss  cone. 
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Figure  III  — 1  illustrates  the  pitch  angle  distribution  of  precipitating  par¬ 
ticles  (at  x=l)  as  the  strength  of  scattering  (the  parameter  ug.)  is  changed.  We 
see  our  expectation  of  increasing  anisotropy  with  increasing  uj  borne  out.  The 
dashed  lines  in  this  Figure  are  plots  of  an  asymptotic  solution  to  Eq .  (6), 
calculated  in  the  limit  that  3f/3x  goes  to  zero: 

fasymp  =  ln(u/ug,)  +  A  (8) 

where  A  is  independent  of  u.  We  see  that  the  asymptotic  form  fits  well  at  pitch 
angles  only  a  little  outside  the  loss  cone,  implying  that  spatial  derivatives 
are  already  small  at  that  point. 

The  variation  of  the  pitch-angle  distribution  along  the  field  line  is 
illustrated  in  Figure  1II-2 .  This  shows  the  filling  in  of  the  loss  cone,  which 
starts  nearly  empty  at  one  end  of  the  field  line  (t=0),  as  t  increases  and  the 
opposite  end  of  the  field  line  is  approached.  We  have  assumed  reflection  sym¬ 
metry  about  the  equatorial  plane,  so  that  the  point  t=1/2  corresponds  to  the 
equator.  This  allows  us  to  construct  the  populations  in  both  loss  cones  at  one 
point  in  space,  by  putting  back-to-back  the  two  distributions  calculated  at 
complementary  values  of  t:  t  and  1-x .  This  has  been  done  in  the  insert  in 
Figure  I1I-2 ,  for  the  case  x=3/4.  This  shows  a  nearly  empty  upflowing  (v||CO) 
loss  cone  along  with  a  fuller  downflowing  ( v h >0 )  one,  much  like  the  observed 
pitch-angle  distributions^ .  None  of  this  structure  could  have  been  calculated 
using  the  earlier  bounce-averaged  theory. 

III. 1.3  Anisotropy  in  the  Continuous  Aurora 

In  the  theory  as  it  is  formulated,  the  diffusion  coefficient  remains  a  free 
parameter.  If  it  were  possible  to  specify  the  amount  of  turbulence  as  a  func¬ 
tion  of  position,  frequency  and  wavenumber,  it  would  be  possible  to  calculate 
the  diffusion  coefficient  using  quasilinear  theory^.  In  practice,  such 
detailed  information  about  the  turbulence  is  never  available.  Instead  we  can 
turn  the  problem  around:  from  the  detailed  observations  of  particle  anisotropy, 
we  can  Infer  the  diffusion  coefficient,  and  obtain  valuable  information  con¬ 
cerning  the  turbulence.  This  empirical  method  has  seen  wide  use  in  the 
radiation  belts  starting  with^.  An  application  of  the  method  to  data  from  a 
post  break-up  aurora,  coupled  with  quantlinenr  eut  I  m.-i  t  «■«;  of  t  },<■  diffusion  <oef 
ficient,  allowed  Lyons ^  to  identify  the  turbulence  there  to  be  1 1 >< ■  result  ol 
unstable  electron  cyclotron  harmonic  waves.  We  go  tin  to  apply  the  method  to 
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I i  t  I  i  I  <■  nliitt'r  7nt  I  >iii<!  In  I  In-  i'oii  (  1  mif  mri  .in  r<>  r;t  f  r  om  flu*  I S  f  S— 7  K.'il  f*1  1  1  to,  kindly 
provided  by  James  Sharber  (Florida  institute  oi  Tecluiology.)  . 

The  technique  is  to  take  the  pitch-angle  distribution  measured  at  one 
energy  during  one  spin  of  the  satellite  and  fit  our  theoretical  pitch-angle 
distribution  to  it  by  adjusting  the  parameter  uj,  •  The  results  of  this  fitting 
are  illustrated  in  Figure  III-3.  This  shows  the  TSIS-2  data  at  several 
energies,  for  one  spin  cycle  at  inv.  lat.  A  =  68.6°.  We  see  that  anisotropy  is 
small  at  low  energies,  E  4.1  keV,  but  that  it  grows  as  particle  energy 
increases.  Using  appropriate  values  of  the  loss  cone  opening  angle  and  bounce 
period,  we  then  determine  4D>b  from  u£  using  Eq.  (7).  The  results  for  the 
pitch-angle  diffusion  coefficient  are  shown  in  Figure  III-4.  Here  we  plot  *D>b 
at  several  energies  as  a  function  of  inv.  latitude  across  the  continuous  aurora. 
We  find  that  falls  with  energy  roughly  as  a  power  law,  E-n,  with  n  between 

one-half  and  unity.  At  fixed  energy,  the  diffusion  coefficient  peaks  at  inv. 
lat.  A  =  69.6°,  where  the  precipitating  energy  flux  into  the  ionosphere  is  also 
largest J.  Cyclotron  emissions  are  also  known  to  be  correlated  with  fluxes  of 
1-10  keV  electrons*’  in  this  way,  suggesting  that  such  waves  may  be  responsible 
for  the  scattering’’. 

1 1 1 . 1 . 4  Discussion 

We  have  demonstrated  how  the  effect  of  particle  loss  out  of  a  loss  cone  in 
velocity  space  is  gradually  smoothed  out  along  the  field  line  by  pitch-angle 
scattering.  Our  technique  was  applied  to  determine  the  amount  of  scattering 
necessary  to  account  for  the  observed  anisotropy  of  electrons  precipitating  in 
the  continuous  aurora. 

The  diffusion  of  electrons  by  turbulent  fields,  such  as  we  have  described, 
is  only  one  aspect  of  the  auroral  precipitation  problem.  Another  side  is  the 
generation  of  the  fields  as  the  result  of  instability  of  features  in  the  par¬ 
ticle  distribution'’.  Each  approach  should  form  part  of  a  self-consistent 
description  of  the  total  system  of  both  particles  and  waves  -  a  challenging 
problem  for  future  work. 
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III. 2.  Ion  Acceleration  in  the  Suprauroral  Region 

It  is  becoming  more  widely  accepted  that  the  energetic  ion  conics^ 
observed  below  shock  structures  in  the  suprauroral  region  are  produced  as  the 
result  of  ion  acceleration  by  the  VLF  turbulence  observed  there^®.  The  tur¬ 
bulence  is  generated  through  the  instability  of  the  auroral  electron  distribu¬ 
tion  accelerated  parallel  to  the  geomagnetic  field  by  the  shock^.  A  model  for 
the  formation  of  ion  conics  in  this  way  was  proposed  by  Chang  and  Coppi^O. 
Acceleration  nearly  perpendicular  to  the  field  line  by  VLF  turbulence  near  the 
lower  hybrid  frequency  is  followed  by  the  adiabatic  folding  of  velocities  as  the 
ions  mirror  and  travel  up  the  geomagnetic  field  line,  creating  the  conic  velo¬ 
city  distribution.  Detailed  calculations  of  conics  using  a  Monte  Carlo  tech¬ 
nique  to  model  the  wave-particle  interaction  were  carried  out  by  Retterer  et 


1)  A  Monte  Carlo  Model  for  Ion  Acceleration 

In  the  lower  suprauroral  region  the  electron  cyclotron  frequency  is 
generally  larger  than  the  electron  plasma  frequency.  Such  an  electron-ion 
plasma  can  support  lower  hybrid  waves^2  with  frequencies  near  the  ion  plasma 
frequency.  For  these  modes,  ^i^^LH^^e*  k||<Mc^,  and  vte/^et<.k-^t.<.vt;:j/fij[ .  It 
has  been  shown  that  these  lower  hybrid  waves  can  be  generated  quite  efficiently 
by  the  electron  beams  produced  by  field-aligned  DC  potential  drops  during  magne¬ 
tic  substorms^,20,23,  Because  of  the  broadband  nature  of  the  turbulence^  >  we 
cannot  use  a  theory  of  acceleration  by  a  coherent  wave^.  Instead,  we  apply  the 
quasi-linear  diffusion  formulation  for  wave  particle  interaction.  Because  of 
the  frequencies  and  wavelengths  involved,  we  may,  in  a  first  approximation,  use 
the  unmagnetized  expression  for  the  diffusion  coefficient  for  the  ions,  given 
for  example  by  Eq .  (10.23)  of  Ichimaru^S.  (The  excitation  of  LH  waves  with  the 
unmagnetized  approximation  dropped  has  been  studied  by  Basu  et  al.26.)  The 
dynamical  picture  is  completed  by  including  the  effects  of  the  geomagnetic  field 
and  the  DC  electric  field,  along  with  the  quasi-linear  term,  in  a  kinetic 
equation  for  the  slow  evolution  of  the  ions: 
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In  this  equation  f  is  the  ion  distribution  function,  v n  Is  the  velocity  along  s, 
which  is  a  coordinate  denoting  position  along  the  field  line,  and  E  and  B  are 
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the  static  fields.  (Drifts  perpendicular  to  B  are  Ignored.)  The  quasilinear 
diffusion  coefficient  D  is 

D  =  (£)2  f—  f-^K  7=7?  SE(k,(D)  TT6(o>-k.v)  (10) 

^  n  J  2ir  J  (2tt  )3  |k|2 

where  Sp;  is  the  spectral  energy  density  in  the  turbulence. 

Because  of  the  Inherent  importance  of  its  spatial  structure,  the  time- 
dependent  solution  of  the  kinetic  equation  for  the  ions  would  be  fully  a  four- 
dimensional  problem  -  too  complicated  to  be  solved  by  standard  finite-difference 
techniques.  Instead  we  adopt  a  particle  simulation  model,  in  which  the 
stochastic  effects  of  wave-particle  interactions  are  described  using  a  Monte 
Carlo  technique.  Because  the  number  of  accelerated  ions  is  small,  we  treat  them 
as  'test  particles'  in  externally  imposed  fields,  instead  of  calculating  the 
fields  sel f-cons Istently . 

From  an  Initial  distribution  in  velocity  and  space,  the  calcuation  of  the 
evolution  of  the  distribution  proceeds  by  following  the  motion  of  a  large 
number  of  ions  with  time.  Because  we  are  interested  in  changes  occurring  on 
scales  much  larger  than  the  size  of  a  gyroorbit,  we  need  not  integrate  the 
equations  of  motion  in  detail.  Instead,  we  follow  only  the  motion  of  the  par¬ 
ticle  guiding  centers  along  the  field  line.  Between  the  velocity  perturbations 
caused  by  interaction  with  the  waves,  it  is  assumed  that  the  ions  travel  in  the 
inhomogeneous  geomagnetic  and  DC  electric  fields  with  constant  energy  and  first 
adiabatic  invariant. 

The  wave-particle  interactions  are  taken  into  account  in  the  following  way. 
In  each  time  step  At,  the  velocity  of  each  particle  is  perturbed  by  an  increment 
Av  chosen  according  to  a  probability  distribution  P^t(Av).  In  kinetic-theory 

terms,  the  resulting  change  in  the  distribution  f  is  given  by  the  Smoluchowski 

">  7 

equat ion- ' 

f ( v^» t+A t)  =  / f(v-Av_,t)  PAt(Av^)  d^A  v_  .  (11) 

To  model  the  quasilinear  diffusion  process,  one  chooses  PAt(^v)  to  be  a  gaussian 
function  with  t he  following,  moments 

Av  )  =  Vv  '  D  At  and  (.  A  vAv^  =  2  J5  At  (12) 
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where  is  the  quasilinear  diffusion  tensor.  With  this  choice,  the  Fokker- 
Planck  expansion  of  Eq.  (11)  yields  the  quasilinear  diffusion  equation  in  the 
limit  that  At+0.  In  the  calculation,  the  probability  distribution  P/^t(Av)  is 
sampled  using  random  numbers.  Among  many  other  applications,  this  Monte  Carlo 
technique  has  been  used  to  introduce  collisions  into  plasma  simulations3^’  39 . 
We  go  on  now  to  apply  the  technique  to  the  ion  conic  problem. 


III. 2.1  The  Formation  of  Ion  Conics 

To  study  ion  conics  in  the  lower  suprauroral  region,  we  carry  out  a  simu¬ 
lation  using  the  following  parameters.  Over  the  altitude  range  from  1000  to 
5000  km  we  follow  the  evolution  of  IT*"  and  0+  ions,  whose  initial  density  distri¬ 
bution  decreases  with  the  second  power  of  the  altitude  -  an  approximation  to 
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Maeda's  model  of  ionospheric  density.  Their  initial  velocity  distribution  is 
an  isothermal  Maxwellian  distribution,  with  a  temperature  of  1  eV.  As  the  simu¬ 
lation  progresses,  we  allow  a  steady  state  to  be  established  by  replacing  every 
particle  which  leaves  the  simulation  by  a  particle  picked  at  random  from  the 
primordial  distribution.  For  these  first  calculations,  we  include  no  DC  poten¬ 
tial  drops. 

The  last  parameter  we  must  specify  is  the  diffusion  coefficient.  To 
evaluate  the  quasi-linear  term,  we  need  to  know  the  spectral  density,  Sg,  of  the 
lower  hybrid  turbulence.  A  self-consistent  calculation  of  Sjr(k,w)  is  compli¬ 
cated  by  the  difficulty  of  determining  the  saturation  mechanism  of  the  lower 
hybrid  instability,  whether  it  occurs  by  nonlinear  evolution  or  simply  by  con¬ 
vection  out  of  the  beam.  We  can  sidestep  that  problem  by  using  the  experimen¬ 
tally  observed  amplitudes  of  the  lower  hybrid  waves,  which  can  range  up  to  50 
31 

mV/m  and  beyond  .  In  order  of  magnitude  we  have,  in  the  range  of  resonant 
velocities, 
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where  |ew |  is  the  observed  amplitude  of  the  wave.  The  resonant  range  extends 
from  a  few  times  the  ion  thermal  velocity  up  to  a  velocity  vxmax  “  u^e  ( "k n / k x ) 
max^H  *  Ufoe  (me/m-[)^^,  where  u^e  is  the  velocity  of  the  electron  beam.  Thus 
ions  should  stay  in  resonance  until  they  reach  energies  near  the  electron  beam 
energy,  1-10  keV.  Above  the  upper  limit  v^max,  should  behave  asymptotically 
as  vi"3.  Unfortunately,  lack  of  knowledge  of  the  wavenumber  spectrum  of  the 
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turbulence  prevents  us  from  specifying  the  velocity  dependence  of  the  diffusion 
coefficient  in  any  more  detail  than  as  a  constant  within  the  resonance  limits. 

In  practice,  we  will  find  that  saturation  of  the  heating  is  not  caused  by  the 
resonance  limit  but  instead  by  convection  out  of  the  region  containing  the 
turbulence.  We  consider  several  cases  for  the  altitudinal  range  of  turbulence: 

from  narrowly  spread,  over  a  few  km,  to  widely  spread,  over  a  thousand  km.  The 
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actual  spatial  distribution  is  poorly  known,  although  evidence  shows  that 
broadband  lower  hybrid  turbulence  will  be  found  below  the  electrostatic  shock 
region  whenever  electron  beams  are  detected. 

Let  us  consider  the  model  with  a  narrow  range  of  turbulence  first.  In  this 
case  the  heating  and  folding  processes  work  almost  independently  of  each  other. 

In  the  steady  state  we  find  that  the  resonant  ions  have  been  heated  to  a  charac¬ 
teristic  energy  of  10-20  eV  when  the  turbulence  range  is  10  km.  Their  pitch- 
angle  distribution  remains  sharply  collimated  at  all  altitudes,  as  it  folds  from 
near  90°  at  1000  km  to  about  150°  at  5000  km. 

In  contrast,  when  diffusion  occurs  over  a  wide  range  of  altitude,  the 
resulting  conic  structures  are  more  diffuse.  The  continued  transverse  heating 
destroys  the  collimation.  Figure  III-5  illustrates  the  ion  conic  formed  in  this 
case,  with  a  scatter  point  plot  of  particle  kinetic  energies  parallel  and  per¬ 
pendicular  to  the  magnetic  field.  The  lower  panel  of  the  figure  presents  the 

particles  In  the  altitude  range  from  1000  km  to  2000  km  -  the  LH  acceleration 
region  in  the  simulation;  the  upper  panel  gives  the  ion  conic  after  it  has  left 
the  acceleration  region,  using  the  particles  from  2000  km  to  3000  km.  The  acce¬ 
lerated  ions  create  a  hot  tail  on  the  ion  energy  distribution^ as  Figure  III-O  she 
Because  of  the  wide  extent  of  the  acceleration  region,  particles  reach  high 

energies  -  up  to  a  few  keV  in  this  case.  The  increase  in  energy  is  not  greater 

because  the  heating  is  sel f-limi t ing :  as  a  particle  gains  energy,  it  moves  out 
of  the  heating  region  more  quickly.  The  broader  conics  calculated  in  this 
second  model  are  much  like  the  conics  seen  by  S3-3^. 

II J. 2.3  Discussion 

We  have  illustrated  the  process  of  conic  formation  through  transverse 
heating  by  lower  hybrid  waves  and  propagation  along  the  geomagnetic  field. 
Although  we  explicitly  considered  only  H+  ions  in  the  calculations  presented 
here,  other  calculations  show  that  0+  ions  can  he  heated  to  energies  comparable 
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to  those  of  H+;  the  number  of  resonant  0+  ions  in  a  hydrogen  dominated  plasma, 
however,  will  be  small. 

Focusing  on  the  broadband  nature  of  LH  turbulence,  we  have  here  adopted  the 
quasilinear  description  of  the  acceleration  process.  Other  workers  have  con¬ 
sidered  stochastic  heating  in  the  presence  of  a  single  wave.  In  the  lower 
suprauroral  region,  where  turbulence  is  intense  but  broadband,  it  is  clear  that 
neither  approach  is  completely  satisfactory.  What  is  needed  is  a  theory  of 
strong  turbulence.  Once  the  wave-particle  interaction  in  such  a  theory  can  be 
described  in  a  probabalistic  way,  it  is  amenable  to  the  Monte  Carlo  treatment 
presented  here. 

Considerable  uncertainty  remains  in  the  model,  however  because  of  the  dif¬ 
ficulty  in  estimating  the  rate  of  the  wave-particle  interaction  process. 
Empirical  estimates  of  the  velocity  diffusion  tensor  based  on  observed  wave 
amplitudes  suffer  because  the  lack  of  wavenumber  measurements  prevents  us  from 
determining  the  phase  velocities  of  the  waves.  We  cannot  rely  on  linear  calcu¬ 
lations  to  give  us  the  wave  spectrum  either,  because  there  appears  to  be  a  real 
difficulty  in  linearly  exciting  waves  of  small  enough  phase  velocity  so  that  the 
resonant  portion  of  the  ambient  ion  distribution  can  account  for  the  observed 
number  of  particles  in  the  conics.  In  addition,  the  self-consistent  evolution 
of  the  wave  spectrum  has  been  ignored  in  previous  work. 

2)  Plasma  Simulation 

To  address  these  problems,  a  plasma  simulation  was  performed  to  provide  an 
independent,  self-consistent  means  of  studying  the  generation  of  the  turbulence 
and  the  resulting  ion  acceleration.  The  suprauroral  situation  was  modeled  by 
allowing  a  weak  (nb/n0  -  10-2),  energetic  (Eg  =  1  keV),  warm  electron  beam 
traveling  along  the  magnetic  field  to  destabilize  a  cool  electron-ion  plasma  (Te 
=  Tj  =  1  eV).  We  set  the  direction  of  propagation  of  the  waves  to  be  nearly 
perpendicular  to  the  magnetic  field,  with  cos^Bg  =  me/tnj  to  reflect  the  observed 
wave  spectral  peak  near  l.r>  times  the  lower  hybrid  resonance  frequency.  The 
velocity  of  the  1  keV  electron  beam  projected  onto  the  propagation  direction  is 
then  about  32  times  the  initial  Ion  thermal  velocity.  The  phase  velocities  of 
waves  excited  by  this  beam  will  be  far  out  on  the  tail  of  the  ion  velocity 
distribution,  where  few  ions  can  resonantly  interact  with  them. 
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Nevertheless,  a  finite  fraction  of  the  ions  are  significantly  accelerated 
in  the  course  of  the  simulation.  We  found  that  tails  of  energetic  ions  formed, 
emerging  from  both  sides  of  the  initial  distribution  at  about  three  times  the 
ion  thermal  velocity;  some  ions  are  accelerated  to  velocities  comparable  to 
those  of  the  electron  beam.  In  addition  to  the  tails  the  core  of  the  velocity 
distribution  showed  evidence  of  nonresonant  beating.  It  can  be  fitted  by  a 
Maxwellian  velocity  distribution,  in  which  changes  in  the  fitted  thermal  velo¬ 
city  reflect  the  changes  in  the  total  wave  energy.  But  the  tails  account  for 
most  of  the  energy  transferred  to  the  ions  in  the  course  of  the  instability: 
immediately  following  wave  saturation,  at  t  *  400  they  already  contain  3% 

of  the  ions  and  account  for  half  of  the  ion  energy. 

The  interpretation  of  these  results  is  clear,  because  it  follows  from  the 

extensive  work  devoted  to  the  high-frequency  analogue  of  the  problem:  electron 
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tail  formation  in  strong  Langmuir  turbulence  .  The  intense  VLF  waves  linearly 
excited  by  the  beam  parametrically  decay  into  lower  phase  velocity  VLF  waves  by 
coupling  through  nonresonant  quasimodes  which  are  driven  to  finite  amplitude  in 
the  turbulent  state.  These  lower  phase  velocity  VLF  waves  are  then  Landau 
damped  by  the  plasma,  accelerating  the  ions  perpendicular  to  the  magnetic  field 
and  the  electrons  (because  of  their  restricted  perpendicular  mobility)  parallel 
to  the  field. 

Several  calculations  support  this  interpretation  of  our  simulation.  First 
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an  analysis  of  the  nonlinear  dispersion  relation  for  the  coupling  of  two  lower 
hybrid  waves  through  nonresonant  quasimodes  was  performed.  Using  the  amplitude 
and  other  parameters  of  one  of  the  linearly  excited  waves  in  our  simulation  as  a 
pump  wave,  we  calculated  the  phase  velocities  of  the  sideband  waves  excited  by 
the  three-wave  coupling  process.  We  found  that  thes ?  velocities  agreed  well 
with  the  velocities  at  the  points  where  the  tails  emerge  from  the  background 
distribution,  supporting  the  argument  that  Landau  damping  of  these  sidebands 
accelerates  the  tons.  Second,  direct  simulation  of  the  parametric  decay  of  a 
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lower  hybrid  pump  wave3  was  found  to  lead  to  perpendicular  ion  acceleration 
accompanied  by  parallel  electron  acceleration.  Finally,  we  formulated  a  simple 
set  of  kinetic  equations  containing  the  quasilinear  equations  with  mode  coupling 
terms  to  describe  the  parametric  processes.  Numerical  solution  of  these 
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equations  produces  ion  velocity  distribution  with  high  energy  tails  and  a  wave 
spectrum  similar  to  the  ones  observed  in  the  simulation. 

III. 2. 3  Conclusion 

Scaling  of  our  simulation  results  to  suprauroral  conditions  gives  results 
which  agree  well  with  data  from  observed  ion  conics:  fraction  of  accelerated 
ions  - 10“3  to  10“2;  average  energy  -50  eV:  maximum  energy  -1  keV:  observed 
lower  hybrid  wave  amplitude  ranging  up  to  50  mV/m.  We  conclude  that  the  VLF 
turbulence  generated  below  field  aligned  potential  drops  in  the  suprauroral 
region  can  account  for  the  acceleration  of  ions  observed  in  ion  conic  events. 
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IV .  PROBE  THEORY 


IV. 1.  Introduction 

Probe  theory  has  been  an  active  area  of  research  for  almost  two  decades 
starting  with  the  original  study  of  Mott-Smith  and  Langmuir^.  Earlier  studies 
dealt  mostly  with  the  laboratory  plasmas.  With  the  advent  of  the  space  age, 
the  theory  had  to  be  extended  so  that  it  can  be  applied  to  space  plasmas  as 
well.  An  interesting  topic  of  space-related  probe  research  is  the  spacecraft 
charging  due  to  emission  of  positively  charged  particles.  We  are  particularly 
interested  in  the  experimental  study  by  Cohen,  Sherman  and  Mullen^ .  Their 
results  on  the  variation  of  the  extended  probe-to-payload  potential  difference 
(4>t)  with  respect  to  the  electron  density  and  to  the  beam  current  (IQ)  need  to 
be  explained  analytically. 

Since  their  observations  indicated  that  the  vehicle  potential  was  indepen¬ 
dent  of  the  neutral  density  and  the  vehicle  pitch  angle  the  effects  of  colli¬ 
sions  and  of  the  magnetic  field  seem  to  be  unimportant.  Thus,  as  a  first 
attempt,  it  is  reasonable  to  consider  the  collisionless,  non-magnetic  field 
electrostatic  theories^.^, 5  to  seek  an  explanation  of  the  observed  results. 
However,  two  important  differences  from  these  theories  can  be  noted  immediately 
(i)  the  existence  of  the  beam  current  and  (ii)  the  different  geometry  of  this 
problem.  We  can  show  that  the  observed  saturation  value  of  <j>t  can  be  explained 
in  terms  of  fairly  simple  considerations. 

The  main  difference  that  is  expected  to  be  due  to  beam  emission  is  the 
limitation  on  the  minimum  value  the  probe  potential  can  attain.  (Positive  ion 
emission  makes  the  probe  attain  a  negative  potential  $p  with  respect  to  the 
ambient  plasma).  In  the  experiment^,  the  beam  emission  energy  is  ~2  keV.  If  <t>p 
were  to  go  below  (-2)  kV,  the  beam  particles  would  not  be  able  to  escape,  and 
this  constitutes  an  obvious  limit  on  $p,  |<f>p|  2  kV.  As  the  beam  current  is 

increased,  |<j>p|  would  keep  increasing  until  it  reaches  this  saturation  value. 
The  extended  probe  is  at  a  distance  of  152  cm  (taken  perpendicular  to  the 
cylindrical  axis  of  the  probe).  If  the  effective  sheath  were  smaller  than  this 
distance,  the  extended  probe  potential  would  reflect  the  ambient  plasma 

potential  (=0),  and  ♦  t“*frl“$p  would  be  the  same  as  (-$p)  =  |<t>pj-  Under  these 
conditions,  <J>t  would  also  saturate  at  2  kV.  The  observed  value  of  ~1  kV 
suggests  that  only  part  of  the  potential  drop  occurs  at  a  distance  ri  such  as 
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that  of  the  extended  probe,  and  the  sheath  size  may  be  larger  than  this 
distance.  A  simple  estimate  for  the  sheath  size  r0  is  obtained  from  the  "free 

2 

fall"  condition  I0  =  r0n0vte  where  n0  is  the  ambient  charged  particle  density 

and  vt  is  the  thermal  velocity  of  the  ions  in  the  ambient  plasma  (beyond  the 
sheath).  For  the  parameters  of  this  problem  r0»ri,  and  the  potential  variation 
from  the  payload  to  the  extended  probe  is  essentially  governed  by  the  Laplace's 
equation  since  the  charge  term  is  negligible  near  the  payload. 


If  we  consider  the  rocket  to  be  an  idealized  sphere,  with  radius  rp, 


and  the  ratio 


P 


H  =  .  _  Ip 
lip i  ri 


This  ratio  is  independent  of  the  beam  current  or  the  saturation  potential.  For 
this  experiment,  depending  on  how  one  defines  the  equivalent  sphere  radius,  rp, 
we  find  P^O.!  to  0.75.  This  provides  a  saturation  value  for  t0  1*5  kV. 


If  we  consider  the  rocket  to  be  a  long  cylinder,  and  use  the  two- 
dimensional  picture  in  its  neighborhood, 

4> p  ~  B£n(rp/r0)  ,  <f>i  ~  B*n(ri/r0)  , 

and  the  ratio 


_  ^n(rL/rp) 

P  *n(r0/rp) 

Now  the  ratio  depends  mildly  on  the  beam  current  through  rQ.  For  the  para¬ 
meters  of  this  experiment,  one  finds  p“0.5  and  the  saturation  value  for  if’t”!  kV, 
which  is  in  good  agreement  with  the  observed  value. 


IV. 2.  Detailed  Theory 

It  is  rather  remarkable  that  the  above  simplified  picture,  without  any 
detailed  matching  of  the  near  and  far  solutions  and  without  using  the  charge 
term,  gave  such  a  good  agreement  with  the  experimental  result.  It  is  clearly 
important  then  to  develop  a  more  detailed  theory  of  this  phenomenon  along  these 
1  i  n  e  s  . 
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In  order  to  develop  a  theory  for  the  mixed  geometry,  it  is  helpful  first  to 
obtain  a  detailed  understanding  of  the  purely  spherical  and  purely  cylindrical 
situations.  The  first  refers  to  a  spherical  probe  and  a  spherical  sheath,  the 
second  refers  to  an  infinite  cylindrical  probe  surrounded  by  a  cylindrical 
sheath  surface.  While  the  basic  differential  equation  for  both  cases  and  some 
analysis  and  numerical  solutions  in  each  case  are  already  available  in  the  work 
of  Lam^ ,  we  have  gone  much  beyond  that  effort  and  succeeded  in  obtaining  very 
accurate  analytical  representations  for  the  space  potential.  These  results  also 
allow  us  to  develop  a  theory  for  the  mixed  geometry. 


where  F  is  a  scaled,  dimensionless  potential  and  x  =  rp/r  where  r  is  the 
distance  from  the  center  of  the  probe  and  rj  is  slightly  larger  than  r0,  the 
free  fall  sheath  size  given  by  I0  =  irr2n0vte,  where  n0  is  the  ambient  charged 

particle  density,  vt  the  thermal  velocity  of  the  ions  in  the  ambient  plasma  and 
I0  the  beam  current. 

Lam3  showed  that 


4 


lim 

T+l 


F(t)  =  (j) 


4 

(t-1)3  , 


(2) 


lim 

T  >» 


F(x) 


1.9x 


2.7  , 


(3) 


and  also  provided  a  table  of  numerical  values  for  F(x)  for  x=l  to  x=10,  and  a 
graph  of  F(x)  over  a  smaller  range  of  x. 

Introducing  x=ex,  x=logx  simplifies  Eq.(l)  to 
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where  prime  denotes  differentiation  with  respect  to  x. 


Obtaining  a  series  solution  for  small  x  is  straightforward;  one  finds, 


F 


+  b^x  +  b2x^ 


(5) 


ao  =  (3/2)*/3  ,  bl  =  |  ,  b2  =  ^  ,  b3  =  g~ . 

Such  a  solution  was  also  obtained  recently  by  Leadon  et  al^,  except  that  their 
coefficients  were  in  the  form  A^  =  a0,  A2  =  acbi,  A3  *  a0b2,  etc.  and  were 
expressed  numerically  in  decimal  form.  Thus,  the  simple  rational  character  of 
b^  and  b2  was  not  recognized.  We  further  noted  that  there  were  simple  relations 
between  successive  coefficients. 

3  3  3 

b2  =  10  bl  ’  b3  K  12  b2  ,  H  -  b3  .  (6) 

which  lead  to  a  simplified,  explicit  equation  for  F, 

F'  -  F  =  g(x) 

g(x)  »  ~  aQ  xl/3  ^1  -  (7) 

An  explicit  solution  in  terms  of  incomplete  gamma  functions  is  thus  obtained, 

F(x)  -  j  a0^YA/3(x)  -  ^  Y 7/3(x)J  .  (8) 

The  series  solution  (5)  is  not  valid  for  large  x,  whereas  (8)  is  well  defined 
for  all  x.  It  was  also  found  to  be  numerically  accurate  to  within  8  parts  per 
million  for  x^l  (when  compared  to  the  exact,  numerically  obtained  F) .  It  can  be 

shown  also  that  for  large  x,  the  result  is  about  7-%  too  low  and  it  reduces  to 

4 

the  form  (3),  with  slightly  different  coefficients.  In  spite  of  all  the  success 
and  convenience  of  the  solution  (8)  it  is  still  only  a  (good)  approximation  and 
one  can  seek  iterative  procedures  that  provide  improvements. 


Some  exact  results  can  be  derived  from  (4),  which  by  two  integrations,  and 
using  the  boundary  conditions  F(o)  =  0,  F^Co)  =  0  can  be  converted  to  the 
integral 

r  x  ,  y.'  l 

F  =  ex  /  e~x  dx'  f  dx"  -  .  (9) 

Jo  /F(x") 

This  leads  to  the  exact  result,  for  large  x, 

F  =  Aex  +  B  ,  (10) 

where 

,  oo 

A  =  /  dx  e-x  F-1/2(x)  ,  (11) 

j o 

f  00 

B  =  -  /  dx  F-1/2(x)  .  (12) 

•'  o 

An  interesting  procedure  to  evaluate  A  is  to  use  the  small-x  expansion  (5)  for  F 
in  (11);  this  leads  to  A  ”  1.913  by  keeping  the  first  few  terms  in  F.  Thus  we 
have  shown  that  the  small-x  expansion  and  the  large-x  behavior  of  F  are  con¬ 
nected  in  a  precise  fashion.  This  is  a  fairly  general  result  for  the  Poisson 
equation  when  the  charge  density  term  is  a  function  of  the  potential. 

An  alternate  approach  to  obtain  the  large  x  behavior  of  F  is  to  assume  a 
direct  series  expansion  in  inverse  power  of  t  to  solve  (1).  Such  a  procedure 
was  carried  out  by  Leadon  et  al^;  this  method,  however,  cannot  show  the  connec¬ 
tion  between  the  small-x  and  large-x  behavior,  as  through  (11)  or  (12).  In 
fact,  an  auxiliary  matching  of  the  two  series  solutions  in  some  intermediate 
domain  of  x  becomes  necessary,  and  the  accuracy  of  determination  of  the  coef¬ 
ficients  of  the  large-x  series  depends  on  how  the  matching  is  carried  out.  We 
use  the  large-x  series  solution,  and  determine  its  two  arbitrary  constants  A  and 
B  by  exploiting  (11)  and  (12). 

It  should  be  noted  that  (8)  can  be  made  more  accurate  by  keeping  higher 
order  terms  of  g(x)  in  (7),  and  if  we  do  so,  the  improved  (8)  precisely  repre¬ 
sents  the  large-x  behavior  with  the  proper  coefficients  A  and  B. 

In  summary  we  have  developed  several  new  results  for  the  spherical 
geometry,  and  a  very  accurate  explicit  eolution  has  been  provided  for  the  poten- 
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tial .  The  methods  developed  are  quite  general  and  can  be  carried  over  to 
cylindrical  geometry  as  shown  below. 


where  G  is  a  dimensionless,  scaled  potential.  Lam  provided  a  table  and  a  graph 
for  G,  by  numerical  integration  of  (13)  over  a  small  range  of  t.  Interestingly, 
as  t+1,  G(t)  has  the  same  behavior  as  F  in  Eq .  (2).  Lam  did  not  provide  any 
analog  of  Eq .  (3) . 

We  have  analyzed  (13)  by  methods  similar  to  the  spherical  case,  and  we  find 
in  terms  of  x=£nT, 


(14) 


G(x)  =  a0x4/3  (1  +  bjx  +  b2X^  +  ...)  ,  (small  x)  ,  (15) 

ao  =  (3/2 )4/3  ,  bl  =  _  ,  b2  =  . 

G(x)  =  Ax  +  B  +  ...  (large  x)  ,  (16) 

A  =  f  dx  e~x  G_1/2  ,  (17) 

'  o 

f  00 

B  =  -  dx  x  e“x  G'1/2  .  (18) 

Jo 

Results  (17)  and  (18)  are  exact,  and  if  we  use  (15)  to  represent  G  in  those 

equations,  we  find  A  «  2.09  and  B  “  -0.75  from  the  first  few  terms  in  (15)  .  We 

have  also  solved  (13)  directly  in  terms  of  a  large  x  series  of  the  form 
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CO 


oo 


3m 

2 


n+1 


(19) 


G(x)  =  AX  +  B  +  Z 


m=l  n=o 


-mn 


a-mx 


where  A  and  B  are  already  known  from  (17)  and  (18)  and  the  constants  cran  are 
expressible  In  terms  of  A  and  B. 

The  accuracy  of  our  results  for  the  cylindrical  case  is  comparable  to  that 
of  the  spherical  case.  Since  Lam  gave  the  numerical  (exact)  results  for  G  for 
only  a  small  range  of  t,  we  cannot  ascertain  this  directly  over  a  wider  range. 
However,  we  have  matched  the  small-x  and  large-x  series  and  found  that  they 
agree  very  well  in  the  intermediate  domain.  An  analog  of  (8)  for  the  potential 
in  the  cylindrical  case  can  also  be  written  down,  but  more  terms  are  required 
for  accuracy. 


IV . 2 . 3  Mixed  Geometry 

The  more  realistic  problem  of  the  space  probes  requires  us  to  consider  a 
cylindrical  picture  when  r«X./2,2.  being  the  length  of  the  cylinder,  a  spherical 
picture  when  r»£/2,  and  the  problem  of  matching.  The  simplest  approach  is  to 
consider  the  plane  of  symmetry  passing  through  the  center  of  the  cylinder  and 
perpendicular  to  its  axis.  In  this  plane,  one  can  use  (19)  for  the  small-r 
(large-x)  solution  and  (5)  or  (8)  for  the  large-r  (small  x)  solution.  This  is 
a  reasonable  approximation,  since  the  role  of  the  charge  term  in  the  Poisson 
equation  becomes  quite  weakened  near  the  probe  cylinder  (i.e.  for  large  t),  and 
only  the  first  two  terms  in  (19)  representing  a  Coulomb  potential  and  a  fixed 
bias  term  are  really  important.  Further  improvements  are  obtained  by  keeping 
the  next  few  terms  in  (19).  While  the  structure  of  the  solution  is  not  changed, 
the  numerical  values  of  A  and  B  are  no  longer  as  in  (17)  or  (18).  These  have  to 
be  determined  now  by  matching  (19)  and  (5)  over  the  intermediate  domain  at  two 
points . 


IV. 3.  Comparison  with  Experiments 

The  saturation  value  of  $t*  in  Figure  5  of  Cohen  et  al^  is  well  explained 
by  our  approach.  The  development  of  the  theory  for  the  mixed  geometry  allows  us 
to  obtain  results  for  intermediate  values  of  t,  which  provides  a  comparison  with 
Figure  A  of  their  work. 
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It  Is  clear  from  the  geometry  of  the  experimental  set-up  that  the  potential 
drop  quite  near  the  long  cylindrical  payload  will  be  better  described  by  the 
two-dimensional  Laplace's  equation.  This  seems  to  be  borne  out  by  the  better 
agreement  of  p  (2-dira)  with  the  observed  value.  Another  test  of  this  picture 
would  be  to  examine  the  variation  of  p  with  current  Ic.  Since  rQ  increases  as 
1/2 

IQ  ,  p  (and  <j>t)  will  decrease  with  increasing  current  beyond  saturation.  This 
can  be  easily  tested  if  additional  experimental  data  become  available. 

To  understand  the  variaton  of  <J>t  for  the  full  range  of  IQ,  we  replotted 

1/4 

that  data  on  a  log-log  plot.  We  find  that  increases  only  as  ~1Q  in  the 

low-current  domain.  This  is  to  be  contrasted  with  for  a  Lam-type  theory^. 

The  two  curves  (observed  data  and  Lam-type  theory)  cross  at  Io~50pA,  <t>t~700V. 

One  must  explain  why  the  observed  curve  remains  above  the  Lam-theory  for  low 
currents  and  why  it  remains  below  the  Lam-theory  for  high  currents. 

Since  Lam-theory  does  not  include  beam  emission,  it  does  not  lead  to  a 
saturation  for  <}) p ,  and  the  increasing  difference  from  the  observed  results  for 
larger  currents  is  due  to  this  limitation  of  that  theory.  A  modified  theory, 
which  takes  Into  consideration  the  beam  emission  and  the  effect  of  the  ion  space 
charge  in  the  emitted  beam  is  required.  We  have  calculated  the  amount  of  posi¬ 
tive  charge  Q+  in  the  beam  up  to  rQ  and  find  it  becomes  *  comparable  to  Q_,  the 
total  negative  charge  on  the  rocket,  for  high  currents.  Working  out  the 
electrostatics  for  this  "tadpole"  charge  distribution  might  provide  a  smooth  and 
gradual  approach  to  saturation. 

Another  physical  effect  missing  in  Lam-type  theories  is  the  "trapped  ion” 
effect.  It  has  been  shown  that  when  the  probe  size  is  small  compared  to  the 
sheath  size,  ion-trapping  will  occur  and  this  will  reduce  the  pull  exercised  by 
the  negative  charge  of  the  rocket.  This  in  turn  will  allow  a  build-up  of  a 
larger  probe  potential  for  the  same  beam  current.  This  effect  might  explain  why 
the  observed  potential  is  larger  than  the  Lam-theory  in  the  region  of  lower 
currents . 
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